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1 Abstract: A spectral decomposition method is used to obtain solutions to a class of 

nonlinear differential equations. We extend this approach to the analysis of the fractional 
form of these equations and demonstrate the method by applying it to the fractional Riccati 
equation, the fractional logistic equation and a fractional cubic equation. The solutions 
reduce to those of the ordinary nonlinear differential equations, when the order of the 
fractional derivative is a = 1. The exact analytic solutions to the fractional nonlinear 
differential equations are not known, so we evaluate how well the derived solutions satisfy the 
corresponding fractional dynamic equations. In the three cases we find a small, apparently 
generic, systematic error that we are not able to fully interpret. 

Keywords: fractional calculus; nonlinear fractional differential equations; spectral 
1 decomposition; operators; eigenvalues 


12 1 . Introduction 

13 Herein we propose a spectral method for solving fractional nonlinear rate equations of a certain kind. 

14 The method is not perturbative, but neither is it exact, since it gives rise to systematic deviations of 

15 the analytic solution from the numerical solution at intermediate times that reaches a maximum value 

16 of 2%. On the one hand, the spectral method provides a remarkable good approximation to numerical 

17 calculation. On the other hand, the source of the small but systematic deviation from the numerical 

18 solution remains a mystery. This paper presents the approach in detail and introduces a new problem 

19 that requires explanation. 
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Despite the advances made into the understanding of complex nonlinear systems in the last half of the 
twentieth century, many physical phenomena failed to be described using the tools of ordinary calculus. 
Nonlocal distributed effects and memory effects observed in relaxation phenomena [1], living systems [2, 
3], wave propagation in porous materials [4] have been more successfully modeled within the framework 
of the fractional calculus [5]. Fractional differential equations (FDE) have been adopted to explain these 
and other complex phenomena [6,7]. Since exact solutions to the majority of FDEs are not available, 
the search for appropriate analytical and numerical methods is a subject of ongoing research. Recently, 
a number of approaches devoted to solving EDEs have been proposed. Examples include Adomian 
decomposition method [8], homotopy perturbation method [9,10], the fractional sub-equation method, 
and the Haar wavelet method [11], to name but a few. However, the convergence region of solutions 
obtained with these algorithms is rather small. 

It was hypothesized [12] that the spectral decomposition method can be extended to the analysis 
of a class of nonlinear fractional differential equations (NEDEs). Herein we demonstrate the method 
by applying it to the fractional Riccati equation (ERE), the fractional logistic equation (EEE) and a 
fractional cubic equation (ECE), where the fractional-order is in the range 0 < a < 1. The solutions 
obtained are shown to have the correct short-time and long-time behaviors. Additionally, they reduce to 
the well known solutions of the ordinary nonlinear differential equations, when the order of the fractional 
derivative is a = 1. In the cases considered herein the exact analytic solution to the ENDE was not 
known previously, we evaluate how well the derived solutions satisfy the corresponding NEDE using 
numerical techniques. In all cases we find a very small, but systematic deviation of the analytic from the 
numerical solutions that has eluded our best efforts to interpret. One possibility, of course, is that the 
numerical technique used to solve the NEDE is the culprit, since it was based on numerically solving 
linear fractional equations. But this remains to be investigated. 

In Section 2 we introduce the spectral decomposition of the solution to define the eigenvalue problem 
for integer-order linear and nonlinear rate equations, as well as NEDEs. In Section 3 we obtain a 
series expansion over the spectrum of eigenvalues and eigenfunctions for the solution to three NEDEs, 
where the exponentials in the solutions to the integer-order equations, also obtained, are replaced with 
Mittag-Eeffler functions (MEEs). Exact solutions to NEDEs are rare in the literature [ 13,1 4] , so to test the 
validity of the analytic results we numerically evaluate how the solutions obtained satisfy the appropriate 
NEDE. To our surprise the error function measuring this fit is not zero, but varies in time, increasing as 
at early times and decreasing as at late times, and reaching a maximum difference of less than a 
few percent at an intermediate time. This non-monotonic scaling difference is shown to occur with the 
solutions to the ERE, the EEE, as well as, the ECE, all with the same qualitative behavior in the error. 
In Section 8 we draw some tentative conclusions including the speculation that this systematic deviation 
may be generic. 

2. Spectral decomposition 


2.1. Integer operator 
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Let us begin by establishing the nomenclature used in the study of the nonlinear differential equations. 
Consider the one-dimensional first-order differential equation 

|.Y(() = OA'(«), (1) 

where X(t) is the dynamic variable of interest and is a generic operator acting on X (t) . Allowing Eq. 
(1) to describe any dynamical system of interest entails the formal solution 

X{t) = (2) 


where xq = X(0) defines the initial condition in the phase space for the dynamic variable and the 
operator Og acts on the initial condition. The exponential operator is formally defined by the series 
expansion 


Oot ^ (^ 0 ^)^ 

&r(fc + i) 


(3) 


so that the solution Eq.(2) can be expressed as 


x(t) = E 


(Opt)" 




«r(;= + i) 




(4) 


where the operator Og^ acts solely on the initial condition. Note that for a linear equation with a constant 
rate A the operator is given by 


Og% 



k 


Xg = X^Xg, 


(5) 


which when inserted into Eq.(4) and summing the series yields the exponential solution to the scalar rate 
equation 

X{t) = e^^xg. ( 6 ) 

It is apparent that Eq.(5) has a form suggestive of an eigenvalue equation and that the solution to the 
general integer-operator rate equation can be expressed as an eigenfunction expansion over the spectrum 
of eigenvalues 

OO 

(t) = ^Ck(j)k (xg) Xk (t) . (7) 

k=0 

The quantity 4>k (a:o) Xk {t) is the eigenfunction, factored into a piece determined by the spectrum 
of eigenvalues {A^; fc = 0,1, 2, ■■}, a piece determined by the initial condition xq, and the expansion 
coefficient Ck determined by the dynamics and overall initial normalization. Inserting Eq.(7) into (1), 
allows us to separate out the time-dependence of the components of the expansion 

j^Xk it) = XkXk it) Xk it) = ( 8 ) 

Correspondingly, the eigenvalue equations are given by 


Og(pk (a^o) — Afc0fc (xo) 


(9) 
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57 and the eigenvalues are determined by the form of the operator. 

In the linear case just considered the operator is the same as before, so the equation for the 
eigenfunction is 

, d(j)k , , 

UXq 

with the solution 

yk 

(i)k (Xo) = . 

The linear eigenvalue spectrum is degenerate Afc = A and the coefficients are determined from the initial 
condition to satisfy 

OO 

=1. 

k=0 


58 The resulting solution is, of course, given by Eq.(6). 


59 2.2. Non-integer (fractional) operator 

Now assume that this general form of a solution to a differential equation translates to the fractional 
calculus domain. Thus, we replace Eq.(l) with the fractional differential equation 


where 0 < a < 1. We assume the fractional derivative to be defined in the Caputo sense: 

1 /■* X'(t) 




a 


) Jo (t- ry 


7 dr. 


( 10 ) 


( 11 ) 


where X'{t) denotes the derivative of X{t) with respect to its argument. Eq.(lO) can be solved 
analytically in terms of the MEE by employing the spectral decomposition introduced above in which 
case we have for the components of the eigenfunction 


fffiXkit) = XkXk {t), 

to obtain the MEE evaluated over the spectrum of eigenvalues 

xk{t) = K {Xkn . 

The MEE is defined by the series [15,16] 


U = 


E 

k=0 


r (ka + 1) 


Consequently, inserting the MEE into the expansion for the solution yields [12] 


OO 

X{t) = 'Y^Ckfk (xo) (Afcf° 

k=0 


( 12 ) 


(13) 


(14) 


(15) 


60 and the eigenvalues are determined by the operator in Eq.(9). The MEE reduces to an exponential 

61 function when a = 1, reducing the series expansion to the ordinary solution of Eq.(2) in that case. 

62 Note that we can adopt the same formal spectral decomposition employed for the integer-derivative 

63 case discussed above for the fractional-order dynamics considered here. But, before we explore the 

64 fractional case, let us examine an integer-order nonlinear dynamic equation. 
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65 3. Riccati equation 


We illustrate the spectral decomposition method, applied to a nonlinear rate equation, using the 
reduced form of the Riccati equation [17]: 

j^X{t) = l-X\t) = OX{t). (16) 

The Riccati equation is put into the form of Eq.(l) by introducing the phase space operator O : 

(H) 

and the formal solution to the Riccati equation can be written as the eigenmode expansion 


X{t) = 


(18) 


k=0 


66 The dependence of the solution on the initial condition can be described by the eigenfunctions equation 

d(t)k (xo) 


(2:0) = [1 - 2:0] 
= \k4>k (a^o) 


dx 


0 


(19) 


The solution to Eq.(19) is given by 


ln0fc (xo) = yin 


1 + Xq 

1 - Xo 


67 where we need to determine the spectrum of eigenvalues {A^} . 
This last expression is inserted into Eq. (18) to obtain 


x{t) = Y.Ct 


fc =0 


1 + Xq 
1 - Xo 




( 20 ) 


It is straightforward to obtain the eigenvalues Xk = —2k by inserting Eq.(20) into Eq.(16) and equating 
coefficients of time dependent terms, as well as the coefficients for equal k values 

C*! 


The series solution is then 


Ck = 2{^] ;k>0. 


X{t) = 2j2i^y 


k=0 


1 + Xo 


The final expansion coefficient is determined to be Ci = —2 from the initial condition, so that the 
solution to the nonlinear initial value problem becomes 


X{t) = 2Y, 


k=0 


Xo + 1 


( 21 ) 
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Additionally, for xq = 0, Eq.(21) reduees to the well-known solution of Rieeati equation [17] 


X{t) = tanh(f). 


For initial values xq 7 ^ 0 the solution to Eq. (16) is obtained by summing the series in Eq.(21) to obtain 


_ tanh(t) + Xq 
Xq tanh(t) + 1 


( 22 ) 


also in agreement with the solution obtain by Davis [17]. 

Sinee we have an exaet solution we ean test the series expansion eomputationally. The error is defined 
as the differenee between the time derivative of the series solution obtained using Eq.(21) and the known 
exaet solution given by Eq.(22) inserted into the right hand side of the Rieeati equation: 


A(t) = 2^ 

fc =0 


/ Xq - 1 
\Xo + 1 


k 

{-2k) - 1 + 


tanh(f) + Xq 
Xq tanh(f) + 1 


(23) 


The sum is estimated by its first 100 terms yielding a difference on the order of 10 which matches 
machine precision, affirming that the series given by Eq. (21) is a valid solution to the Rieeati equation. 


4. Fractional Rieeati equation 


4.1. Spectral decomposition solution 


We now establish that a generalization of the spectral decomposition method leads to valid solutions 
for certain NFDEs. We start from the fractional form of the Rieeati equation 

Xx(t) = 1 - X^t) = OXX) (24) 

where 0 < o; < 1. If the fractional derivative is of the Caputo type, Eq.(24) has the operator given by 
Eq.(17). Consequently, the solution to the ERE is of the form 

X{t) = sum’^^QCk(l>k (xo) Ea (Afcf), (25) 


and the eigenvalue problem is the same as for the integer-order Rieeati equation, that is, given by Eq.(19). 
Consequently, we have the same spectrum of eigenvalues obtained in the integer-order RE, obtained at 
early times using the stretched form of the MEF. In the same way the expansion coefficients are the same 
as previously obtained in order to satisfy the initial condition and we obtain for the solution to the ERE: 


X(t) = 2j2 

k=0 


Xq 


Xo + 1 


Eo, {-2kE) - 1 . 


(26) 


We note that the solution obtained to the ERE does not coincide with that obtained by Zhang et al. 
[14]. However we point out that these authors use Jumarie’s modified form of the Reimann-Eiouville 
fractional derivative [18], which explicitly satisfies the Eeibniz condition: 




[/WsWI = sM^ 


+ /(^) 


d'"g{t) 
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D“X(t)=1-X"(t) 



D"X(t)=X(t)-X"(t) 



D“X(t)=-X(t)-X"(t) 



10^ lO' 10'* 10’ 10' 

Time, t 


Figure 1. Solutions to the fractional differential equations, (left) Fractional Riccati 
differential equation, as defined by Eq. 26. (middle) Fractional logistic equation, as defined 
by Eq. 37. The growth rate is A = 1.00. (right) The cubic fractional equation, as defined by 
Eq. 48, with a = b = 1.00. The order of the fractional derivative in all cases is a = 0.75. 
Continuous lines correspond to solutions obtained with the spectral decomposition method, 
dashed lines correspond to numerical integration of given equations. The integration step is 
h = 10-3. 


for the derivative of the product of two analytic functions f(t) and g(t) . It is evident that the Caputo 
fractional derivative does not satisfy the Eeibniz condition [15,16] and consequently, although the 
starting dynamic equations look the same, that because of the restrictions on the fractional derivatives, 
the problem solved by Zhang et al. [14] is different from the ERE considered here. 

We see in this figure that the analytic solution does extremely well in following the numerical solution 
to the ERE, but it is not exact. The numerical solutions for various initial conditions, obtained using 
Eq.(26), where the infinite sum is approximated by its first 100 elements, are demonstrated in the first 
panel of the tryptic in Figure 1. One can easily verify that the proposed solution satisfies the initial 
condition and has the correct long time behavior. But since there are no known exact solutions to the 
initial value problem for the NRE, we use the numerical test introduced previously to test the veracity of 
Eq.(26). 


4.2. Testing the solution 

Since the exact analytic solution to the ERE Eq. (24) is not known, one possible approach to check the 
validity of Eq.(26) is to check that this solution satisfies the fractional differential equation in question. 
Applying the Caputo definition of fractional derivative to the analytic solution to the RME yields 

LHS = —X(t)^2Y^I-^^) (-2k)EA-2kn, (27) 
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Figure 2. The differenee A(t) between the RHS and LHS of the fraetional Rieeati equation, 
when the solution is assumed to be Eq. 26. Conseeutive panels eorrespond to an inereasing 
order of the fractional derivative a. Color lines correspond to range of initial conditions, as 
denoted by the legend in lower right plot. 


while the right hand side of the ERE is 

xq-I 

Xq + I 

Again the potential error is defined in terms of the difference 

A(f) = RHS - LHS, 


(-2kH) - 1 


RHS = = 1 - [2^ 


k=0 


(28) 


(29) 


which for an exact solution should be zero. This error variable is estimated numerically, and its 
behavior is plotted in Eigure 2 for four values of a. It is apparent that the difference variable A(f) 
departs from zero at early times, gradually increasing, reaching a maximum value at t Ri 1, and finally 
decreases at long times. The behavior of A(f) and its maximum value Amax < a few percent, clearly 
depends on the order of the fractional derivative, a, and on the initial condition xq. The precise form of 
power-law scaling of A(f), present at short and long times, is obtained adopting an approximation to the 
MEE. At early times, the series defining the MEE (Eq. 14) reduces to a stretched exponential function 


=limE„(-Af“) 


1 


r (1 + oi) 


+ ... = exp 


XH 

r(l + a)J ’ 


(30) 
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Figure 3. Short and long time power law behavior of the error function A(f). Both regimes 
are approximated by equations 32 and 33, respectively. The order of the fractional derivative 
a = 0.90. 


while at late times the MLF has inverse-power law behavior 


= lim {-XE) = 


t—^OO 


Ar(l - a)' 


(31) 


Thus, at early times the difference A(f) scales as and after some algebra this difference is determined 
to be 


2a 


f 

A°(t)=liinA(t) = 

t^o 1 (1 + a) 


(1 


(32) 


while the long-time behavior of the difference is determined by direct calculation to scale as t 

^ ''xq + 1' 


t 


A-(f) = lim A(f) = 

t^oo 1 — a) 


2 log 


loa-2 ( ;eo±1) 

I ™ I 1 V 2 J^-a 

-h Xo -h i-=-7-- 

T [1 — a) 


(33) 


Figure 3 illustrates the validity of both approximations to the difference A(t) for selected orders of the 
fractional derivative. It is evident that the deviation of the solution obtained by spectral decomposition 
deviates from the exact solution in systematic ways at both late and early times. The maximum difference 
for a = 0.9 is shown in this figure to be below 2%, but the explanation as to why the deviation behaves 
the way it does remains elusive. Therefore we examine the solutions obtained by applying the technique 
to two other well known nonlinear initial value problems extended to the fractional domain. 


5. Fractional logistic equation 


5.1. Spectral decomposition solution 

Next we examine the fractional logistic equation (FLE) 

Aa'(() = A“X(() (1 - A'(f)) 


(34) 
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where X{t) is the normalized population and A is the population growth rate. As before we define a 
phase space operator 

O = (1 - x) ^ (35) 

ox 

and write the formal solution as Eq.(25). The eigenvalue equation is slightly different from that for the 
FRE 

(xo) = A"xo (1 - Xo) 

dxQ 

= ^k(t)k ( 2 : 0 ) (36) 


After some algebra we obtain for the eigenfunetions 

0fc(^o)= " 

VI - ^0/ 

and after some algebra, using the speetral deeomposition method, we arrive at the eigenvalues \k = 
— fcA", the expansion eoeffieients Ck = Cf and from the initial eondition Ci = —1. Inserting these 
quantities into the formal solution Eq.(25) yields 

= E^{-k\-r). (37) 


where again xq is the initial eondition.The numerieal solutions for various initial eonditions, obtained 
using Eq.(37), where the infinite sum is approximated by its first 100 elements, are demonstrated in the 
middle panel of the tryptie in Figure 1 . As you might expeet the solutions are qualitatively the same as 
in the first panel. 

When a = 1, Eq.(37) reduees to the well-known solution to the logistie equation 


X{t) 


__ 

Xo + (1 -xo)e-^*’ 


(38) 


yielding sigmoidal population growth from the initial value X (0 ) = xq to the saturation level X (oo) = 1. 

Note that even though the FRE and the FEE have quadratie nonlinearities, the speetra in the two ease 
are different. The former inereasing as 2k and the latter as k, whieh also determines the differenee in the 
normalization parameter Ci in the two oases. 

The solution, given by Eq.(37), is identioal to the one obtained by one of the authors [13] using the 
Carleman embedding teehnique. The Carleman embedding approaeh has sueeessfully determined the 
solution to many nonlinear differential equations as presented and disoussed by Kowalski and Steeb 
[19]. The teehnique was applied to a FNDE in [13] to obtain an infinite-order hierarohy of fraotional 
moment equations, using Eaplaoe transforms and solved using a matrix diagonalization method. A 
modest deviation of the analytie solution from the numerieal integration of the FEE was noted in this 
earlier analysis. 
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Figure 4. The differenee A(t) between the RHS and LHS of the fraetional logistie equation, 
when the solution is assumed to be Eq. 37. Columns eorrespond to deereasing values of of. 
a = 0.90 (left), a = 0.75 (middle) and a = 0.50 (right). The initial value Xq = 0.75 in 
all eases. Top row demonstrates the error funetion A(t) for a range of growth rate A values. 
Middle row demonstrates the effect of rescaling time, while bottom row achieves complete 
overlap of A(t) curves generated for different A through rescaling ?/—axis variable. 

110 5.2. Testing the solution 


Here again we check if the solution given by Eq.(37) satisfies the EEE by examining the difference 
equation. The left hand side of Eq.(34) is 


LHS 

111 while the right hand side is 


dt°‘ 


A'(i) = E 

k=0 


Xq - 
Xo 


1 


k 

(-fcA") K (-A:A“t"), 


(39) 


RHS 


A“X(f)(l-X(f)) 




k=0 


3^0 - 1 
Xo 


(-kre) [i-Yl 


k’=0 


3 ^ 0-1 

3^0 


k' 




(40) 

(41) 


112 The difference or potential error is defined as before by Eq.(29). 

113 Eigure 4 demonstrates the behavior of A(t) for selected values of the order of the fractional derivative 

114 a, the initial condition xq and the unperturbed growth rate A. We observe the same dependence of the 
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differenee A(t) on time as in the case of the FRE. Note again that the the difference is not monotonic, 
but has a maximum value of less than 1%, depending on a and xq. The location of maximum in time is 
a function of A. Rescaling both the x—axis and y—axis as follows 


t —y tX 

A{t) A(f)A-“ 


(42) 

(43) 


normalizes the A{t) curves, making them independent of the growth rate, as shown by the bottom row 
of curves in the figure. 

Adopting the short-time and long-time approximations to the MLF we find that at early times the 
difference A{t) scales as : 

j-2a. 

A°(f) = JimA(f) = (^0 - 1)' , (44) 


while the long time behavior scales as f “ 


A“(f) = lim A{t) 

t—^oo 



(xo - 1 - logxo) 


log^xo 

r(l-a) 


(45) 


The accuracy of the algebraically determined scaling is numerically demonstrated in Figure 5. Here 
again, the deviation of the solution determined using the spectral decomposition technique to solve the 
FLE from the exact numerical solution, is non-monotonic and never exceeds 1%. Thus, the worst case 
using this technique is better than the best obtained using many approximation techniques. 


6. Cubic fractional differential equation 


6.1. Spectral decomposition solution 


The final equation that we consider is a fractional differential equation with a cubic term (FCE) 

^X{t) = -aXit) - hX\t) = OXit). (46) 

The formal solution given by the spectral decomposition requires the solution to the eigenvalue equation 

(a^o) = -Xq [a + 6xo] 

uXq 

= Xk(j)k (xo) ■ (47) 


The algebra providing the solution to the eigenvalue problem will be presented elsewhere. Here we 
record that the solution is given by sum over the eigenvalue spectrum 


x(‘) = E 

k=0 


(2k-1)!! 
(2k)\\ 


b 2 

a-^O 


Xo 


^0 + 1 


--E, 


-xl + l 


(2k -f 1) at°‘), 


(48) 


where we have adopted the notation of double factorial (2k)\\ = 2k (2k — 2) (2k — 4) ■ ■ ■. The solutions 
obtained using Eq. (48) are presented on the right panel of Fig. 1. 
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Figure 5. The differenee A(t) between the RHS and LHS of the FLE. Conseeutive panels 
eorrespond to an inereasing order of the fraetional derivative a. Color lines eorrespond to a 
range of initial eonditions, as denoted by the legend in lower right plot. The growth rate in 
all eases is A = 1. 


For the integer-value ease a = 1, the MLF ean again be replaeed by an exponential in Eq.(48) and 
the series summed to yield the exaet integer-value solution 


Xit) 



(49) 


129 6.2. Testing the solution 


Onee again we eheek how well this solution satisfies the NFDE in question by taking the differenee 
between the two sides of the FCE. The differenee or potential error function A(f) is illustrated on Figure 
6 for a = 6 = 1. The short-time behavior of A(f) is 


A°(f) = limA(f) 


t2a 

P (1 + a) 


xl {xo - 1 )^ 



r (1 -f cr) 


(a^o - 1)^ 


while the long-time behavior is 


(50) 


A'^(f) = lim A(f) 

t^oo 


Cl (a,Xo) 


r(l-a) 


+ C 2 (a, Xq) 


4.-20. 


t-3a 


r2(l - a) 


+ C 3 (a, Xq) 


r3 (1 - a) 


(51) 
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Figure 6. The differenee A(t) between the RHS and LHS of the FCE. Conseeutive panels 
eorrespond to an inereasing order of the fraetional derivative a. Color lines eorrespond to 
range of initial conditions, as denoted by the legend in lower right plot. 


where Ci, C 2 and C 3 are coefficients depending on a and xq, and there is no need to record their values 
here. 

Here again the deviation of the two sides of the FCE using the spectral decomposition analytic solution 
and the exact numerical calculation of the fractional derivative results in scaling as depicted in Figure 6 . 
At early times the deviation increases as whereas at late times the deviation decreases as an inverse 
power law in time reaching a maximum at an intermediate time that does not exceed a few percent. 

7. Numerical solutions 

Since the analytic solutions to ordinary nonlinear differential equations are notoriously difficult to 
find, one turns to numerical methods for assistance. Due to the nonlocal property of the fractional 
derivatives, known numerical methods are not easily transferable to the fractional calculus. Herein 
we adopted the Adams-Bashforth-Moulton predictor-corrector technique, developed by Diethelm et al. 
[20,21] to investigate the solutions to fractional Riccati equation, fractional logistic equation and cubic 
fractional equation in purely numerical fashion through the numerical integration of the listed equations 
The main motivation is to compare the numerical solutions with the ones obtained through spectral 
decomposition method. 
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Figure 7. The differenee 5{t) between the numerie integration of fraetional differential 
equations and the solutions obtained with speetral deeomposition method. 


Solutions for selected values of a and xq are presented on Fig. 1, together with the solutions obtained 
with the spectral decomposition. The overlap of both curves in all cases is extremely good. To better 
visualize how close the spectral decomposition solution Xsp{t) is to the numerical solution 
on Figure 7 we plot the difference 


<5(f) = ^NUAlit) — Xspit). (52) 

145 In all cases the difference 5{t) is smaller than 1% and is characterized by an increase in short time scale, 

146 maximum value at intermediate times and and a decrease at large times, where the difference scales as 

147 This result, taken together with the behavior of A(f) demonstrates that the spectral decomposition 

148 method provides analytic solutions that are very close to the true solutions to the fractional nonlinear 

149 differential equations. 

150 8. Discussion 

151 It has been demonstrated that the spectral decomposition solution to fractional nonlinear dynamic 

152 problems with quadratic and cubic nonlinearities systematically deviate from the evaluation of the 

153 fractional derivative in a common way. The systematic deviation is not monotonic, but for both quadratic 

154 and cubic nonlinearity grows as the same power law at early time, decays as the same inverse power law 
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at late times, and has a maximum deviation of the analytic solution from the numerical evaluation of the 
fractional derivative of a few percent at some intermediate time. Only the size of the coefficients change 
with the parameter values and the type of nonlinearity, but not the functional form in time. It appears 
that the scaling form of the deviation is a consequence of the MLF in the spectral decomposition, which 
changes as at early times and as t~°‘ at late times. 

Recall that the FLE was the only other case with which we could compare the solution obtained here 
with that obtained using a different method and they turned out to be identical. However, the hierarchy 
of equations solved using the matrix method to solve the FLE [13], implicitly assumed a version of 
the Leibniz condition. The correspondence of that earlier solution and that obtained using the spectral 
decomposition method suggests that the deviation from the numerical solution, obtained using the latter 
technique is related to the Leibniz condition. Unfortunately we have not been able to identify precisely 
how this comes about and this remains a speculation. 
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